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1. INTRODUCTION 

Later 2014, following the success of European Space Agency (ESA) Copernicus acquisition 
programs with open-source product development, remote sensing processes and the search for new potential 
applications has been an unprecedented success [1]. Among the potentialities offered by this change of 
context, the exploitation of the temporal axis in particular is found where once obtaining a temporal series of 
images on the same site was rare and difficult [2], today many download platforms make it possible to build 
such batteries very quickly and anywhere in the world [3]. In parallel with obtaining these data, solutions for 
treating them with large computing powers are developing [4]. Synthetic aperture radar (SAR) images have 
remarkable temporal properties: as the images do not suffer from cloud cover, and the radar is an active 
system, the signal is stable between two acquisitions under interferometric conditions [5]. Also, among the 
advantages of having long time series, one finds in particular the fact to improve the quality, by simple 
temporal filtering, or by more evolved methods, but also the detection of change or more largely, the 
follow-up of activity. In return, the SAR images, by their coherent nature, are subject to the speckle 
phenomenon [6]. 

The binary detection principle is shown when testing change detection where modified areas are 
looked for between two dates. Even before the change detection decision, having a visualization method of 
the images highlighting this information is useful for an operator [7]. Conventionally, the visualization of a 
change detection scenario involves the use of complementary colored channels (green/magenta, red/cyan, or 
yellow/blue). The two channels encode the respective intensities of the images of the couple. As a result, the 
presence of different behaviors between two image’s results in color, while the pixels of a stable amplitude 
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will appear in gray levels [8], this method can be extrapolated to the case of the only three images, since the 
red, green, and blue (RGB) space is three-dimensional. But for N>3 images, the visualization of the changes 
is not so immediate anymore. It is possible to go through animations, as in [9] which make it possible to 
differentiate the permanent type changes from those which are occasional changes, but these products remain 
video objects [10], which do not address the problem of temporal reduction, allowing the result of a single 
product summarizing all the changes that have occurred [11]. 

In this paper, we seek to find such a representation based on colored composition and that exploits at 
best the time axis for activity monitoring. Several works in the literature propose color compositions that use 
the interferometric coherence information of a temporal stack. Among them, Akbari et al. [12] proposes such 
an approach called multi temporal change (MTC) for the visualization of lava extension zones. Red and 
green channels encode amplitude images of the interferometric pair and the blue channel the level of 
interferometric coherence. This color composition does not really address the reduction of the temporal 
dimension [13]. Simply recommend choosing the pair minimizing both the temporal and geometric baseline. 

In the same way, research [14] shows several RGB color compositions where coded channels are 
either an amplitude or a coherence calculated on a Sentinel-1 interferometric pair at 12, 24 and 48 days of 
time baselines. Other colored compositions simply use as input images of three interferometric coherence. 
Liao et al. [15] proposes a representation called CovAmCoh, again allowing visualizing changes on a pair of 
images under interferometric conditions. Coding of the red channel is doing by using coefficient of variation 
of SAR; here it is the spatial coefficient of variation, and not the temporal one, this parameter then gives 
information on the local heterogeneity [16]. As two images are considered, it is possible to average the two 
locally calculated coefficients of variation maps. The green channel encodes the amplitude, and the blue 
channel encodes the interferometric coherence. All these compositions only illustrate the relevance of the 
information contained in the level of interferometric coherence [17]. 

The first works that really address the use of the temporal dimension are presented in [18]. The 
authors propose, in a first product called level 1a, to separate the intensities related to the dry season coded in 
blue, those linked to the wet season, coded in green. Finally, the red channel encodes the interferometric 
coherence, which makes it possible to bring out in this color the anthropoid targets. Colin-Koeniguer et al. 
[18] propose a RGB visualization product called level 18, “superior” in the sense that it involves in the 
coding a temporal statistical parameter which is that of the variance of intensities. Red encodes this variance, 
the green channel encodes the average intensity, and finally the blue channel is a combination of the average 
interferometric coherence and a coefficient proportional to the intensity. The latter product 18 therefore uses 
the notion of temporal variance. More recently, Office National d’Etudes et Recherches Aérospatiales 
(ONERA) proposed a colored representation of temporal stacks of SAR data [19]. This representation is 
intended to make appear in color the pixels which are the seat of potential changes, and to leave in grayscale 
pixels unchanged. 

This approach differs from previous studies on at least three criteria: the first of which is that a 
crucial temporal parameter used is the temporal coefficient of variation, which has remarkable statistical 
properties. No spatial estimate is made in the composition. The second of which is that the color composition 
passes through the use of the hue-saturation-value (HSV) color representation space before being converted 
into the conventional RGB space. Finally, it only uses, at the input, the amplitude, and can therefore be tested 
on georeferenced data of the Sentinel-1 ground range detected (GRD) type. 

In this article, ambition is twofold. First of all, it is necessary to explain the temporal statistical 
properties of the criterion used in this approach and which partially explain its efficiency; the next step is to 
analyze the results on a global scale, following the implementation of the method on the Google earth engine 
(GEE) platform. This implementation of work is an opportunity to make a return experience of this platform. 
The article is structured, in section 2, the proposed principle of visualization is presented, this work addresses 
the theoretical statistical framework for modeling the coefficient of variation useful for visualization. In 
section 3, these results make it possible both to better understand the effectiveness of the method and manage 
its settings. The limits and potential improvements to the GEE platform are discussed and several types of 
results are analyzed, from which different potential applications result. Finally, in section 4, a conclusion 
allows us to summarize the main results obtained and to discuss the prospects for change. 


2. THE PROPOSED ALGORITHM 

This section, it is explaining the general principle of the method, to switch to the HSV color space, 
where the hue of the H channel is representative of the time, the saturation channel S is used to highlight the 
colors where the change is significant, and the value V corresponds to a notion of classic radar intensity. For 
this reason, our visualization method is called REACTIV for fast variation and detection in radar time-series 
with the use of the factor of variation. In practice, several choices are possible to have color association by a 
specific time agreeing to the change that one wishes to highlight [20]. Indeed, for N dates, several types of 
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temporal events can be considered as changes: the punctual appearance of an object (a vehicle), an 
appearance or permanent disappearance (the construction or the destruction of a building), or more complex 
events such as seasonal variations on plant parcels. Given the diversity of these situations, there is no single 
coding solution, but several solutions to adapt to the types of events to highlight or not [21]. 

Visualization method REACTIV uses intensity (V-channel) the image of maximum intensity: this 
maximizes the chances of seeing the event whatever it is. For hue, a color is associated with the date for 
which the maximum intensity is reached. This choice will be particularly suitable for a punctual event 
(the presence of a boat). In the case of a permanent appearance, the color will be representative of the date 
when the signal will be maximum; it will not necessarily be the date of the appearance of the object, but one 
of the dates at which it is particularly “strong” [22]. 

In the case of permanent destruction, the date will belong to one of the dates when the object was 
visible. It is possible to use in the intensity channel, an average of the images of the stack, in order to benefit 
at the same time of a speckle filtering [23]. However, in this case, all changes of a one-off nature disappear 
rapidly, visualization will above all allow to see the seats of seasonal or persistent phenomena. Finally, a 
saturation could be coded with a factor of temporal variation whose interest is to react to the presence of 
changes. Thus, the saturation of the color will be all the stronger as the temporal signal attached to a pixel 
will have a specific temporal evolution [24]. 


2.1. Dynamic settings 

In order to debug the color coding, the actual parameters used for visualization must be defined 
between 0 and 1 [25]. Here, the choices made on the hue and value parameter are explained. The remaining 
choice, relative to the coefficient of variation, is based on the analysis of the statistical properties of this 
parameter [26]. 

H hue: the colored composition affects a color on a specific date. Built-in batteries are not 
necessarily acquired with regular time sampling [27]: when an acquisition is missing from the stack, it causes 
a jump in acquisition dates. Also, our choice is to set the date linearly with respect to the chosen acquisition 
range [28]. If t/ and f2 are starting, and end dates of the stack, that are going to be displayed, then each image 
acquisition date ¢ is recovered for, in order to calculate the time fraction ft, as shown as in the (1). 


t-t1 


ft= (1) 


~ t2-t1 


This parameter is a real between 0 and 1, it serves directly to code as it is the channel H of the image hue, 
also it should be noted in this composition, that the first and the last date have very similar colors because the 
HSV color palette is continuous and loop on itself [29]. 

The intensity (V): the intensity is coded by the maximum amplitude channel. If there is no change, 
and the temporal periodicity is respected, then this intensity follows a classical law gamma order L [28]. 
Also, the dynamics’ choices made can be done in a conventional way in radar imaging thresholding the 
image on a maximum value of dynamics [30]. For visualization purposes, the values stored in dB were 
converted into linear amplitude levels, then threshold the dynamics to | (0 dB) and applied a contrast factor 
or gamma parameter equal to 1:3. 


2.2. Dynamics of the coefficient of variation: theoretical aspects 

The originality of the REACTIV composition lies mainly in the use of the temporal variation 
coefficient to calculate the saturation parameter. It is this parameter that will or will not make our eye 
interpret the color as related to a change. This section is interested in the statistical properties of this 
coefficient that allow to the better understanding its behavior. 


2.2.1. Coefficient of variation of a speckle area 
In radar images, it is generally believed when the amplitude of a texture less speckle which follow 
the law of a Rayleigh-Nakagami law [31], as the (2): 


JL 2 
= 2 ME Enzi- (7) 

RN[u, LI) = Tp) e (2) 
where u is a shape parameter and L is what is called the number of equivalent looks. Commonly, the law 
used will make estimation spatially on similar zone. This work focuses on the statistics of time. Also, 
assuming that for any pixel which belong to a Rayleigh Nakagami law speckle zone has not been undertaken 
for any variation, the numerous realizations of the amplitudes over a time of a pixel follow the same law. 
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This amount [32] to transforming the hypothesis of spatial periodicity into a hypothesis of temporal 
periodicity. In particular, as in (3): 


_ r5 
Mı = Fira) 


et m = p’ (3) 


This allows to find the coefficient of variation, which is the ratio of the standard deviation and the average as 


in (4): 
o \m2-mi T(L)P(L+1) 
y=is = 2—1 (4) 
u mı T'(L+1/2) 
Equations here demonstrated the truth: For every average, the amplitude of the speckle, factor of the variation 
value could be similar for all the stable zones. So as to estimate this parameter behaves, currently the study in 
the change of the estimator is one of the important goals. Kendall [33] suggest calculating the difference of a 


function g (mı, mz), this can be done by a first order limited extending of the role about mo, and mo, 2. The 
process leads to find this formula, proposed in [32]: 


1 4m3—m3m2+m2m4—4mym2m3 


var(y) = (5) 


4N m4 (m2-m?) 


where N is the amount of temporal pictures which is measured, showing the periods | to 4 of a Nakagami 
law, in this way the change in estimation such as L shown in (6): 


i irot arra-aLr(L+t) -1(L+4)°) 


var(y) = (6) 


4N r(L+ 1/2)*ur(L)2-F(L+4) ) 


In this method where L=1, also \~0. 522723 taken var (X)=0. 137881/N, (A)=0. 3713/VN 
respectively for changes and regular deviancy. Applying a multi-look on a SAR image variation the L 
parameter. Sentinel-1 GRD data is provided with a number of the same looks presented at L=4.9, computed 
for the theoretical swath and mid-orbit value [34]. This value of L gives for the coefficient of variation 
associated with the speckle \=0.2286, var (\)=0. 0216/N and (\)=0.1616/VN. 

The equation presented another superior featured of level: 1, this equation also proved that the 
1 
Ww 
have a constant value through a difference reducing the square root of the number of pictures in the stack. 
Several simulations carried out show that when the, v(’) was computed on the amplitude following to a 
Rayleigh Nakagami law, this lead to return back to the truth of a law of Rayleigh Nakagami. The featured is 
very important to simulate many analyzing results and discussions [35]. 


relationship between v(”) and N is inverse since therefore it demonstrates that factor of variation should 


2.2.2. Coefficient of variation of a speckle profile with temporal rupture 

After having analyzed properties of the coefficient of variation for zones of pure speckle, this work 
proposes to simulate the way in which this coefficient behaves on a speckle zone which has a temporal 
rupture. This simulation is performed: It was set a number N of temporal samples, then generate 10° profiles 
of size N, following a Rayleigh Nakagami law of parameter p=0.3 (-11 dB) and L=4.9. the same number of 
profiles are also generated for which the last value of the profile is replace by fixed amplitude, considering 
scaled magnitudes values between mj+2.5 dB and mı+20 dB, representative of a longer temporal rupture. 
Then these realizations make it possible to calculate a coefficient of variation unique, both for pure speckle 
and for a temporal break [36]. 

This estimate is made for all the N's ranging from 2 to 200 as shown in Figure 1. The value curve 
(A+) also superimposed on curves, for having an idea of the statistical distance margin between the pure 
speckle and temporal rupture. It is found that the coefficient variation of the pure speckle converges in 
function from N to the expected theoretical value (0.22), and for the other targets converges more or less 
rapidly towards values which are a function of the ratio between the amplitude of the target and that of 
speckle. Position of maximum value between the coefficient of variation of population with rupture and that 
without rupture depends on both N and the relative amplitude target/speckle. A remarkable result is that these 
curves do not depend on the average value of the speckle: with a fixed number of images, the coefficient of 
variation, in the presence of a temporal rupture, depends only on the jump of rupture between the speckle and 
the target, it deviates rapidly from the value of the coefficient of variation of a stable speckle [37]. 
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Figure 1. Simulated statistics of the coefficient of variation on a “‘stable”+ speckle type profile + temporary 
temporal rupture 


2.3. The standardization choices 

Locally, the number of images available in the same area varies, to avoid depending on v(’) on the 
number of pictures, it is necessary the normalization distribute with the help of averages [38], and the 
theoretical standard deviations. The statistical law of the coefficient of variation, supposed Nakagami on 
speckle will be approximated by normal law. A normal law does not allow overcoming the dissymmetry 
characteristic of Nakagami law, but to go through this approximation allows us to have a form of first-order 
empirical normalization by the number N [39]. In this context, the following empirical standardization shown 
in the (7) is proposed: 


YEE 
Ve T 0.25 (7) 
A theoretic mean of typical value deviancy of “constant” speckle with L=4.9 as demonstrated in 2.2.1 this 


empirical standardization aims, as shown in Figure 2. A saturation of constant areas is decreased nearby of 
0.25 and spreading the variations to greater saturation. 


Saturation coefficients 


Channels 


Figure 2. Empirical normalization of the saturation coefficient 
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3. RESULTS AND DISCUSSION 
3.1. The results obtained locally and on the platform 

GEE is an online platform for digital terrestrial data processing [40]. Its advantage is to store the 
petabytes of data collected by satellites accessible in open source. Then it makes it possible to implement and 
test treatments through a development interface (JavaScript API) and visualize the results on a global scale, 
benefiting from all the computing power of Google. The REACTIV method has been implemented on the 
platform to be applied to Sentinel-1 data, to demonstrate the value of visualization around the globe, and to 
better understand the limitations of the platform. Our approach consisted in comparing on the same area, the 
results obtained on the platform with those obtained with datasets downloaded and processed locally. 

The study area that has been selected is Saclay in Ile de France. The Sentinel-1 GRD batteries were 
downloaded from the ESA platform. Then the data was pre-processed with sentinel application platform 
(SNAP) software, then refined the image registration with the GeFolki algorithm [41], and finally applied the 
REACTIV method implemented in the MATLAB language. 

The stack consisted of 57 images of descending passes, acquired between 29-01-2016 and 
04-04-2017, with a constant average incidence. Figure 3 shows both the representation obtained on the 
platform, and that obtained locally on downloaded data, for the vertical transmit and horizontal receive (VH) 
polarization. Visually, the representations do not seem to differ. They let appear the parcels cultivated with 
strongly saturated colors, as well as locally, zones of construction of which could be verified that they 
corresponded effectively to building sites set up at this time. 


Figure 3. A zoom on a representation made from the data of the platform and the data downloaded locally, 
then flunked 


3.2. The results obtained: quality registration of statistical dynamics 

Today, on the platform, only georeferenced GRD-type Sentinel-1 data is available. Single look 
complex (SLC) data from Sentinel-1 cannot currently be embedded because earth engine does not support 
images with complex values because of its inability to compute them during the pyramid decomposition of 
images without losing the information of phases. Each GRD scene included in GEE has pretreated with 
SNAP by performing thermal noise suppression, radiometric calibration, and terrain correction performed 
with the shuttle radar topography mission (SRTM) digital elevation model (DEM) 30 or advanced 
spaceborne thermal emission and reflection radiometer (ASTER) DEM for areas greater than 60 degrees of 
latitude. No radiometric correction related to relief is currently provided. 

Then the values of the geometrically corrected images are converted to decibels via logarithmic 
scaling. In order to save memory space, the platform makes the choice to convert float32 values to unsigned 
integers of 2 bytes uint16, and only retains the values between the 1st and the 99th average percentile of the 
values before applying the quantification. This may have the consequence of modifying the original statistics 
of SAR signals. Calibration of the temporal data of the platform and local data verified by the method 
proposed in [42]. It did not bring a substantial gain to the stability of our images, so the radiometric 
calibration based on metadata was kept. We analyzed the statistical distributions of the coefficient of 
variation for different homogeneous areas selected the image. Note that in this case, which allows us to look 
at the distribution of this parameter are spatial samples, even though the parameter has been calculated only 
temporally. The histogram of the estimated coefficient of variation on the local data is shown in Figure 4, for 
an area that appears to belong to a “stable” speckle zone (wooded area). In this case, the distribution seems to 
be well modeled empirically by a law of Rayleigh Nakagami, whose expectation and calculated standard 
deviation. 
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Figure 4. Local statistics of coefficient of variation on a “stable” speckle zone 


By maximum likelihood are respectively equal to X\=0. 21648 and (\)=0. 0214, values close to the 
theoretical values found in section 2.2, respectively 0.2286 and 0.0212. A hypothesis test with Rayleigh 
Nakagami law of these parameters concludes the adequacy of this law. The same estimates conducted on 
GEE give an average of 0.1992 and 0.0200, which provides slightly biased estimates compared to those of 
locally processed data. The statistical analysis was reproduced of three change-speckle zones, selected from 
plots of different colors in Figure 3. 

For these three fields, the average values found for the coefficient of variation change between 0.3 
and 0.4, and the standard deviations between 0.03 and 0.05. According to hypothesis tests put in place, the 
distributions obtained no longer follow the Rayleigh Nakagami law. In summary, locally statistical analysis 
shows that the coefficient of variation in a zone without change follows a Rayleigh Nakagami law with 
theoretical the values found. On the platform, the estimation results are slightly biased, although first-order 
errors are negligible in estimating the coefficient of variation: it is mainly the statistics on permanent 
broadcasters, whose values are in the distribution queue, which will be affected. It should be kept in mind 
that the current choice of a dynamic compression of the platform can affect quantitative results based on the 
statistics, even if qualitatively, visualization is not affected. 

Whether on the platform or locally, the data is georeferenced using SNAP, which ensures their 
matching. However, the accuracy of the result obtained can be further improved. On the locally downloaded 
data, a sub-pixel registration method, which has the necessary accuracies for the interferometric mode, has 
been applied [41]. It leads to finding shifts from one image to another not insignificant, up to 2 pixels, with a 
standard deviation of distance found of 0.4 pixels, these residual offset can have a significant influence on the 
saturation and the hue obtained in the composition, especially around elements high lights: an offset on a 
fixed target may result in a change in the edges of the same target, at the same time, the date of a one-time 
event found can also be changed. Figure 5 show the effects were illustrating on a zoom which made on the 
site of Saclay, where several diffusers surrounded by a colored halo. 


3.3. The results obtained of dynamic applications practices around the world 

While keeping in mind the limitations of the platform, mainly related to the imperfections of 
georeferencing, it allows anyway to analyze very quickly the results obtained anywhere on the planet. 
Although it is intended to highlight the targets appearing punctually. It is seen on the example of Saclay that 
the color composition also allowed to visualize other types of changes, such as agricultural dynamics, and 
construction or destruction of urban elements. These two themes are therefore addressed, before addressing 
the detection of boats. 

One of the first observations concerning color visualization is the marked saturation of colors on all 
agricultural parcels. In fact, this dynamic can give complementary information to traditional optical data. It is 
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therefore likely that this type of visualization adds value to the classification of agricultural parcels. This 
result is generalized to the entire planet explored: all agricultural parcels appear with much higher temporal 
variation coefficients than forest patches, for which seasonal effects are much lower. 

In addition to the presence of cultivated agricultural parcels, the composition colored makes it 
possible to put forward building constructions. This has already been validated on the Saclay site. Another 
example is given in Figure 6, which shows the evolution of a district of Beijing during the year 2017. 

Many colorful spots appearing in the image reflect frequent changes in land use. On a large scale, it 
becomes possible to map all urban areas undergoing rapid development. In the same way, it is conceivable to 
study the traction of structures. In Figure 7, the result of REACTIV representation applied to the city of 
Palmyra in Syria is an example. It shows a large number of modifications, probably related to the destruction 
in 2015 following the capture of the city by the Islamic State. 


Figure 5. Zoom on a representation made on the data of the platform on the left, and data downloaded and 
recalibrated on the right 


28-11-2014 04-02-2016 


Figure 7. The REACTIV composition at the palmyra site, 2015 
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For maritime traffic, the REACTIV visualization is particularly sensitive to targets that appear in a 
timely manner, it allows to identify very quickly all the boats of sufficient size, in all the port and coastal 
areas. In Figure 8, the products obtained around the Singapore port area, for which maritime transport is 
particularly important in the world, are visualized. Also, the temporal visualization makes it possible to put 
forward at the same time, the accumulations of boats on the zones of expectations, but also the maritime 
routes. Given this visualization, it is possible to perform on the ocean s’ zones, temporal detection based on 
the use of the coefficient of variation, with a threshold positioned at the theoretical average value plus the 
standard deviation 0.22+0 1616/VN. The fact of going through a temporal detection only and not space, 
makes it a particularly fast tool. The spatially averaged detection map in Figure 8 was obtained over the 
entire Sentinel-1 acquisition period (from the end of 2014 to the end of 2017) in a few seconds, to obtain a 
map that can be interpreted in terms of density traffic over this period. 


Figure 8. The REACTIV composition in Singapore from 2015 to the end of 2017 on the left, and the density 
map of boats detected on the right in Europe 


4. CONCLUSION 

In this article, a REACTIV color composition, based on a coding of information in the HSV plane 
have been proposed, in order to optimize the analysis of a large SAR time series and in particular to address 
the problem of detection of changes. Temporal coefficient of variation whose statistical properties have been 
studied in the framework of Nakagami law achievements. It has been shown that, a speckle zone, this 
coefficient has a constant mean value which depends only on L, a number of equivalent looks, and that its 
variance decreases with the square of the number of images. A time lapse in speckle profile will result in an 
increase in the average value of this coefficient. The method has been implemented under the GEE platform, 
whose main limitations concern dynamic compression and the perfectible quality of the default registration of 
GRD data, the impact of which will have to be measured. Finally, on a global scale, it has been introduced 
visualization products for different applications that have demonstrated the effectiveness of the method. 
These promising results lead researchers today to consider switching from visualization to the actual 
detection of changes. Another challenge will be to manage the multi-event context, that is to say, to have an 
ability to detect several changes over time in one place. 
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